Gene Co-expression Modules

Bulk vs. Single-Cell RNA-seq

Author

Brian M. Schilder, Bioinformatician II

library(dplyr) 
library(ggplot2) 
root = "./"
source(file.path(root,"MAIN.R")) 
import_parameters(params)

## **** __Utilized Cores__ **** = 2$subsetGenes
## [1] "protein_coding"
## 
## $subsetCells
## [1] 500
## 
## $resolution
## [1] 0.6
## 
## $resultsPath
## [1] "./Results"
## 
## $nCores
## [1] 2
## 
## $perplexity
## [1] 30

load("./Data/monocle3_CDS.RData") 
nCores <- 4

Graph-autocorrelation Analysis

Find genes of interest through graph-autocorrelation analysis.

# pr_graph_test_res = monocle3::graph_test(cds, neighbor_graph="knn", cores=nCores)
# pr_deg_ids = row.names(subset(pr_graph_test_res, q_value < 0.05))

Louvain Modules

• Identify co-expression modules.
  find_gene_modules essentially runs UMAP on the genes (as opposed to the cells) and then groups them into modules using Louvain community analysis.
• Advantages: Much faster than other alogrithms like WGCNA.
  
• Disadvantages: Seems to be limited of ~1000 genes.

gene_module_df = monocle3::find_gene_modules(cds[var.genes[1:1000], ], 
                                             resolution=1e-2, 
                                             reduction_method = "UMAP",
                                             max_components = 3,
                                             cores=nCores, 
                                             random_seed = 2019)
data.table::fwrite(gene_module_df, "./Results/sc.modules_by_gene.txt", sep="\t")
createDT(gene_module_df)

Pheatmap

• Plot the eigengene expression of each gene co-expression module (y-axis) within each cell cluster (x-axis).

# Create pheatmap table
cell_group_df = tibble::tibble(cell=row.names(SummarizedExperiment::colData(cds)),
                               cell_group=monocle3::clusters(cds)[colnames(cds)])
agg_mat = monocle3::aggregate_gene_expression(cds, gene_module_df, cell_group_df)
row.names(agg_mat) = stringr::str_c("Module ", row.names(agg_mat))
colnames(agg_mat) = stringr::str_c("Cluster ", colnames(agg_mat))
data.table::fwrite(data.table::data.table(data.frame(as.matrix(agg_mat)), keep.rownames = "Module"),
                   "./Results/sc.modules_by_cluster.txt", sep="\t" )
# Plot
pheatmap::pheatmap(agg_mat, cluster_rows=TRUE, cluster_cols=TRUE,
                   scale="column", clustering_method="ward.D2",
                   fontsize=6) 

3D Scatter

# plotly R documentation: https://plot.ly/r/
# rgl::plot3d(x = gene_module_df$dim_1, y=gene_module_df$dim_2, z=gene_module_df$dim_3, ) 
p.3d <- plotly::plot_ly(gene_module_df, 
                     x = ~dim_1, y = ~dim_2, z = ~dim_3, 
                     color = ~module,  
                     text = ~paste("Gene:",id),
                     colors = "Dark2",  
                     # alpha = .7, 
                     size = 5
                     ) %>%
  plotly::layout(title = 'scRNA-seq Gene Co-expression Modules')  %>%
  plotly::add_markers() %>%
  plotly::layout(scene = list(xaxis = list(title = 'UMAP 1'),
                     yaxis = list(title = 'UMAP 2'),
                     zaxis = list(title = 'UMAP 3')))
htmltools::tagList(setNames(list(p.3d), NULL)) 

Module Expression

monocle3::plot_cells(cds,
                     genes=gene_module_df,
                     group_cells_by="cluster",
                     color_cells_by="cluster",
                     show_trajectory_graph=F) #  rasterize = T (for super high-res figures)

Bulk RNAseq Modules vs. scRNAseq Modules

• Co-expression modules were identified through the WGCNA analysis of bulk monocyte RNA-seq data, conducte by Katia Lopes.
  
• Here we compare the modules identified in the scRNA-seq data here (sc.modules), to the modules identified in the bulk monocyte data (bulk.modules).

# Prepare bulk modules
bulk.modules <- data.table::fread("./Data/Katia_Modules/Gene_modules.tbl.txt") 
# bulk.modules$ensembl.id <- gsub("\\.*","", bulk.modules$gene_ids)
# ensembl.dict <- ENSG_to_HGNC(ensembl_ids = unique(bulk.modules$gene_ids), 
#                              reference_genome="grch38", 
#                              input = "ensembl_transcript_id")
if(!file.exists("./Data/Katia_Modules/ens.geneid.gencode.v30")){
  system("scp schilb03@chimera.hpc.mssm.edu:/sc/hydra/projects/pd-omics/GRCH38.p12/ens.geneid.gencode.v30 ./Data/Katia_Modules")
} 
ensembl.dt <- data.table::fread("./Data/Katia_Modules/ens.geneid.gencode.v30")
ensembl.dict <- setNames(ensembl.dt$GeneSymbol, ensembl.dt$gene_id)
bulk.modules$gene.symbol <- ensembl.dict[bulk.modules$gene_ids]  
bulk.modules <- bulk.modules %>% dplyr::rename(module=moduleColors)
createDT(bulk.modules)

## Warning in instance$preRenderHook(instance): It seems your data is too
## big for client-side DataTables. You may consider server-side processing:
## https://rstudio.github.io/DT/server.html

# Prepare sc modules
sc.modules <- data.table::fread("./Results/sc.modules_by_gene.txt") %>% dplyr::rename(gene.symbol=id)
createDT(sc.modules)

Module-Module Overlap

• Calculate enrichment scores for all combinations of bulk RNAseq and scRNAseq gene modules.
  
• This tells you whether some modules are detectable in both datasets.

genomeSize <- length(union(unique(bulk.modules$gene.symbol), monocle3::fData(cds)$gene_short_name))

MOD.dt <- module_vs_module_enrichment(mod.df1 = bulk.modules, 
                                      mod.df2 = sc.modules, 
                                      genomeSize = genomeSize, 
                                      verbose = F) 

## [1] "mod.df1 contains 66 unique modules."
## [1] "mod.df2 contains 8 unique modules."
## [1] "Conducting enrichment tests on 528 module-module combinations..."
## [1] "528 enrichment tests conducted in 0.99 seconds."
## [1] "22 enrichment tests were significant (at FDR ≤ 0.05)."
## [1] "19.7% of mod.df1 modules showed enrichment for a mod.df2 module."
## [1] "87.5% of mod.df2 modules showed enrichment for a mod.df1 module."

createDT(MOD.dt)  

Plot Module-Module Overlap

# library(heatmaply)
mod.cast <- reshape2::acast(MOD.dt, 
                            Module1.name~Module2.name, 
                            value.var="Jaccard.index")
hm <- heatmaply::heatmaply(mod.cast, 
                            k_row = 5, k_col = 5, 
                            height = 10, 
                            column_text_angle = 0, 
                            xlab = "scRNA-seq Louvain Modules",
                            ylab = "bulk RNA-seq WGCNA Modules", 
                            cexRow = .7,
                            key.title = "Jaccard Index") # file = "module-module.heatmap.html"

## Registered S3 method overwritten by 'seriation':
##   method         from 
##   reorder.hclust gclus

# shiny::includeHTML("module-module.heatmap.html")
htmltools::tagList(setNames(list(hm), NULL)) 

Module GO Enrichment

• Iteratively test for enrichment of ontological terms in all modules (both bulk and single-cell).
  
• Take a while for long lists of modules (e.g. the 60+ bulk-RNA-seq modules).

mod.enrich.path <- "./Results/gprofiler2.module.enrichment.txt"
if(file.exists(mod.enrich.path)){
  print(paste("Module GO Enrichment file detected. Importing...",mod.enrich.path))
  mod.enrich <- data.table::fread(mod.enrich.path)
} else {
  bulk.mod.enrich <- gprofiler2.module_enrichment(bulk.modules, dataset = "bulk-RNA-seq")
  sc.mod.enrich <- gprofiler2.module_enrichment(sc.modules, dataset = "sc-RNA-seq")
  # Merge
  mod.enrich <- data.table::rbindlist(list(bulk.mod.enrich, sc.mod.enrich))
  # Correct for multiple tests
  mod.enrich <- mod.enrich %>% dplyr::mutate(FDR = p.adjust(p_value, method="fdr"), 
                                             Bonferroni = p.adjust(p_value, method="bonferroni")) %>% data.frame()
  mod.enrich[mod.enrich == "NULL"] <- NA # is.na(mod.enrich) <- mod.enrich == "NULL" 
  data.table::fwrite(mod.enrich, mod.enrich.path, sep="\t") 
}

## [1] "Module GO Enrichment file detected. Importing... ./Results/gprofiler2.module.enrichment.txt"

mod.enrich.sig <- subset(mod.enrich, Bonferroni <= 0.05)
createDT(head(mod.enrich.sig)) 

Module-Module Similarity: Gene-level vs. Term-level

• Compute how similar each of the top conserved module-module pairs are in terms of GO enrichment.

mm.overlap <- module_vs_module_enrichment.terms(MOD.dt = MOD.dt, 
                                                mod.enrich = mod.enrich.sig,  
                                                top.terms = F,
                                                row.limit = F)

## [1] "Comparing overlap of enrichment terms for all module-module combinations."

createDT(mm.overlap)

Histogram

# Prepare dataframe
jaccard.corr <- merge(x = MOD.dt[,c("Module1.name","Module2.name","Jaccard.index")]  %>%
                        dplyr::rename(Jaccard.index.gene=Jaccard.index), 
                      y = mm.overlap[,c("Module1.name","Module2.name","Jaccard.index", 
                                        "Module1.top.term","Module2.top.term")] %>% 
                        dplyr::rename(Jaccard.index.term=Jaccard.index), 
                      by=c("Module1.name","Module2.name"))
# Remove zeros
# jaccard.corr <- subset(jaccard.corr, Jaccard.index.gene > 0 & Jaccard.index.term > 0)
## Histogram
ggplot(jaccard.corr) + 
  geom_histogram(aes(Jaccard.index.gene, fill="Jaccard.index.gene"), alpha=.5) + 
  geom_histogram(aes(Jaccard.index.term, fill="Jaccard.index.term"), alpha=.5)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Scatterplot

• See how well module-module similarity scores (Jaccard index) correlate between the gene-level comparisons vs. the term-level comparisons.
  
• Varying how many of the top enrichment terms to use to calculate term-level module-module similarity (e.g. 10, 50, 100, all) did not seem to signficantly impact the correlation between methods (generaly hovered between Pearson’s r = .4-.5).
  
• Label each pair of modules in the scatterplot using the top ontological enrichment terms for each module respectively.
• Jaccard.sum is the sum of both the gene-level and term-level Jaccard similarity scores.
  • This can be used as a way to identify module-module pairs that were similar at multiple levels (and thus have strong evidence of conservation).

## Correlation scatterplot
jaccard.corr <- jaccard.corr %>%
# label.dat <- subset(jaccard.corr, (Jaccard.index.gene >= .05| Jaccard.index.term >=.2) ) %>% 
  dplyr::mutate(Label=paste0("x: ",round(Jaccard.index.gene,2),
                             "\ny: ",round(Jaccard.index.term,2),
                             "\n\n",
                             # Mod 1
                             "Mod1 name: ",Module1.name,
                             "\nTop Term: ",Module1.top.term,
                             "\n\n", 
                             # Mod 2
                             "Mod2 name: ",Module2.name,
                             "\nTop Term: ",Module2.top.term),
                Jaccard.sum = Jaccard.index.gene + Jaccard.index.term)

corr.stat <- cor(jaccard.corr$Jaccard.index.gene, jaccard.corr$Jaccard.index.term, method="pearson")
# Scatterplot
sp <- ggplot(jaccard.corr, aes(x=Jaccard.index.gene, y=Jaccard.index.term, text=Label, color=Jaccard.sum)) + 
  geom_point() +
  # ggrepel::geom_text_repel(data = label.dat, aes(label=Label)) +
  geom_smooth(method=lm, se=T) + 
  labs(title = "Correlation between gene-level vs.term-level module-module Jaccard similarities",
       subtitle = paste("Pearson's r =",round(corr.stat,3)))
ggp <- plotly::ggplotly(sp, tooltip ="text" )
htmltools::tagList(list(ggp))